Collapse times of dipolar Bose-Einstein condensates 
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We investigate the time taken for global collapse by a dipolar Bose-Einstein condensate. Two semi- 
analytical approaches and exact numerical integration of the mean-field dynamics are considered. 
The semi-analytical approaches are based on a Gaussian ansatz and a Thomas-Fermi solution for the 
shape of the condensate. The regimes of validity for these two approaches are determined, and their 
predictions for the collapse time revealed and compared with numerical simulations. The dipolar 
interactions introduce anisotropy into the collapse dynamics and predominantly lead to collapse in 
the plane perpendicular to the axis of polarization. 
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Long coherence times and a high degree of controlla- 
bility make ultracold atomic gases suited to the study 
of non-equilibrium states of many-particle quantum sys- 
tems. One example is the collapse of a Bose-Einstein 
condensate (BEC) [IllllllllSlElEligby using a Fes- 
hbach resonance to change the s-wave scattering length 
fls from positive to negative. Two limiting cases can be 
identified: global and local collapse depending, respec- 
tively, upon whether the (imaginary) healing length as- 
sociated with as is of the same order as, or much smaller 
than, the size of the BEC. The mechanism of global col- 
lapse is an instability of the monopole collective excita- 
tion mode which grows exponentially, causing the entire 
condensate to collapse in 3D [9]. Meanwhile, for local 
collapse it is high-lying phonon modes whose amplitudes 
grow fastest \\^. Local collapse is expected when there 
is a sudden large change in within a large system [S]. 
The stability of a trapped BEC can be parametrized by 
ks = N as lap m nj ini Hg, where N is the number of 
atoms and Op is the radial harmonic oscillator length of 
the trap. The system collapses when the interactions are 
attractive {ks < 0) and \ks\ exceeds a critical value \kg\. 
An important parameter defining collapse is the collapse 
time, tc, which is the time taken for atomic three-body 
losses to become significant [5 . Meanfield simulations 
including three-body loss have reproduced experimental 
results reasonably well [Ti| . 

The ^^Cr condensates made by the Stuttgart group [TS] 
are the first to have large dipole-dipole interactions (in 
addition to the s-wave). Dipolar interactions are long- 
range and partially attractive, and thus the properties 
of dipolar BECs (DBECs) are rather intriguing. In re- 
cent experiments |TB1 HZ] the collapse of a DBEC was 
triggered by reducing the repulsive s-wave interactions, 
with the collapse proceeding anisotropically and on a 
global scale. Motivated by these experiments, we the- 
oretically model the timescale for the global collapse of a 
DBEC. This is achieved through mean-field simulations 
and two semi-analytic methods. The first semi-analytic 



method is based on plasma physics treatments of the col- 
lapse of electrostatic [HI [H] and electromagnetic ^\ HI] 
wavepackets. We apply it to the case where the initial 
state is weakly interacting. The second method is valid in 
the opposing interaction-dominated Thomas-Fermi (TF) 
regime; essentially we run the usual ballistic expansion 
equations [22", "23^ in reverse. 

The macroscopic wave function describing a BEC 
ip{r,t) satisfies the Gross-Pit aevskii equation (GPE) 



Zm 



dd - 



(1) 



where g — ATrh^as/m and m is the atomic mass. The con- 
fining harmonic potential, Vt{v) — mcj^ (p^ + 7^2;^) /2, 
where — x^+y^ , is cylindrically symmetric, with radial 
(axial) frequency LOp {ojz — l^p)- The term <I>dd(r,t) = 
J d^r'Uddi^ — ^')n{r',t) accounts for the dipolar inter- 
actions Udd{j^) = AioM^(l ~ 3 cos'^ 6')/(47r|r|^), where 9 is 
the angle between the separation vector r — r' and the 
polarization direction, taken to be the z-direction.The 
strength of the dipolar interactions are characterized in 
terms of kd = Nad/ap, where Ud = i-iop?'ni/{12TTfi^) is 
the length scale of the dipolar interactions. The last term 
in the GPE represents three-body loss where is the 
recombination rate. For the semi-analytic methods we 
assume that the DBEC collapses globally, maintaining 
a single-peaked density profile. In the GPE simulations 
some features of local collapse, discussed elsewhere [24], 
are observed. This typically occurs after considerable 
global collapse and as such the dynamics remain in good 
agreement with global collapse predictions. 

Analogous collapse occurs in plasma physics where 
electrostatic [TSJ [TH] and electromagnetic [201 HI] wave 
packets in a turbulent plasma undergo nonlinear self- 
focusing and Zakharov collapse. Following approaches 
for plasma wavepacket dynamics |19j , and work by Lush- 
nikov for DBECs [25], the equation of motion for the 
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mean square radius of the DBEC, (r^) = (p^ + z^), is 
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where (r^) = / r'^\ip\^(Pv and (z^) = / z^j^/ipd^r. We 
assume an initial stable state with total and kinetic en- 
ergy, Ea and Ek^), and, (r^) = (r2)o, (z^) ^ {z^)q with 
d{r'^)o/dt — d{z^)o/dt — 0. Upon changing the strength 
of the s-wave interactions from Oso to the upper bound 
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where E' — Eq — EsQ{l—as / a^Q) . Rearranging, the upper 
bound for the collapse time is. 
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where collapse occurs when (r^)y = 0. 

In the limit of weak interactions the condensate can be 
modeled by a Gaussian ansatz (GA): 



^GA = ^/n/I^I^ 
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The energy E of this ansatz is evaluated using the GPE 
energy functional, and has four contributions: kinetic 
(Ek), trap {Et), s-wave {Es), and dipole (Ed) [26j . 
By minimising E with respect to the radii Ip and / 
PSI E7] the variational solutions for Eq, {r"^) f and 
are found. This enables the time it takes for the BEG to 
go from (r2)o to (r^) j to be evaluated, via Eq. For 
{r'^) / = this defines the collapse time. 

When the interactions become dominant it is then ap- 
propriate to use the TF approximation, wherein the zero- 
point kinetic energy is neglected. In this limit, the non- 
dipolar GPE under harmonic trapping is known to sup- 
port an exact scaling solution given by [l^J [23] , 
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(6) 
(7) 



where no{t) = \hN j^-KRpitfR.Xt)]- This class of solu- 
tion remains exact even in the presence of dipolar interac- 
tions [28', '29' . Substitution into the GPE yields equations 
of motion ([s]) and ^ for the radii [29^. These describe 
the two lowest-lying excitation modes, namely the axis- 
symmetric quadrupole and the monopole excitations. 

It is important to establish how well the TF and GA 
approximations reproduce the exact ground state in the 
parameter space of fc^ and fc^. We define their regimes 
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FIG. 1: The regimes of validity for the GA and TF ground 
states in the parameter space of ks and kd- Above the short 
dashed black curve there are no stable ground states to the 
GPE. The red solid (green long dashed) curve bounds the 
regions of validity of the GA (TF) solution, defined to be 
when their energy is within 5 % of that of the true GPE 
ground state. Considered are the trap ratios of (a) 7^ = 0.1, 
(b) 1 and (c) 10. In each of the figures the circles [squares] are 
the values for fcso and kd that are employed in Figs. 2(a-c) 
and Figs. 3(a-c) [Figs. 3(d-f)]. 

of validity to be when the energy of the GA or TF solu- 
tion differs by less than 5% from the energy of the GPE 
solution. Fig.[l]maps out these regimes for (a) 7^ — 0.1, 
(b) 1 and (c) 10. The short dashed black curve marks 
the threshold for collapse: above it there are no stable 
ground state solutions to the GPE. Hence, for a given kd 
there is a unique fc^ below which the system is unstable 
to collapse. The region bounded by the solid red (long 
dashed green) curve marks the region of validity for the 
GA (TF) solutions. In general we find that the GA gives 
a good approximation to the ground state for weak in- 
teractions and in regions close to the collapse threshold. 
By contrast, the TF solutions are most accurate in the 
presence of strong s-wave interactions and for parameters 
well away from the collapse threshold. Furthermore, trap 
geometry plays a key role in the validity of the solutions. 
In prolate, 7 < 1, traps the GA (TF) approximation has 
a large (small) region in which it is valid, while for oblate, 
7 > 1, traps the opposite is true. This is because a pro- 
late (oblate) dipolar BEG experiences a net attractive 
(repulsive) interaction due to the dominance of end-to- 
end (side-by-side) dipoles, and for such interactions the 
GA (TF) approximation works well. Note that close to 
the collapse threshold zero-point kinetic energy plays a 
significant role in stabilising the condensate and so the 
TF approximation does not provide a good description 
for ground states there. However, the TF approach can 
still be employed to model the collapse dynamics provid- 
ing that the interactions dominate over zero-point kinetic 
energy during the dynamics. In practice this is achieved 
by beginning with an initial state that is well within the 
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FIG. 2; Ap(t) and A2(t) for fc^o = 60, fc^ = and fc^o = kd = 
30, with (a) 7^ = 0.1, (b) 1, and (c) 10. The solid (dashed) 
curves co rresp ond to Ap (A^); thick black curves are solutions 
to Eqs. ( |8|9[ ) and the thin red (medium green) curves are 
the full solutions to the GPE with (without) loss. The red 
triangle indicates when the sudden onset of loss has begun, for 
isTs = 2 X 10"^* cm'^/s [17]. The black square is tc as derived 
from Eq. (j4]| using the GA and the green circle indicates when 
p = for the GPE solution in the absence of losses. For the 
GPE simulations we have defined Ap(t) = {p^(i))/{p^(0)) and 
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regime of stability (a TF initial state), and then sud- 
denly switching to a point in the parameter space that 
is deep in the collapse regime, i.e. bypassing the thresh- 
old for collapse. As the density increases during collapse 
Es + E]j continues to dominate Ek- 

We now study the dynamical collapse in the TF limit. 
We define the scaling parameters Ai(t) = Ri{t)/Rio, 
where Rio are the initial radii (i = p,z). Then, under 
general time-dependent changes in ks(t) and kd(t), the 
time evolution of Xp{t) and Xz{t) is given by, 
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where the unit of time is l/oj^ 



/(«) 



a — 15Koa^/i?po 
The initial aspect ratio. 



kq, is evaluated from Eqs. (8 9) for A^ = [261 129 



Figure [2] shows the evolution of Ap (solid curves) and 
A2 (dashea curves) as calculated from the TF equations 
(thick black curves) and numerical simulation of the GPE 



with (thin red) and without (medium green) three-body 
loss. Specifically, the case where ks is switched from 
fcso = 60 to (fcrf — kdo — 30) is considered, for (a) 
7^ = 0.1, (b) 1 and (c) 10. For these parameters the 
initial state of the BEG spans the regimes of validity of 
the TF and GA approaches. Gomparing the GPE results 
to the TF analysis excellent agreement for the majority 
of the collapse is observed. Significant deviations only 
occur close to the point of collapse (Ap — > 0) when the 
GPE collapse bounces and turns into an expansion of the 
system, consistent with the recent observation of a d- 
wave explosion following collapse [T7|. Importantly, the 
collapse is highly anisotropic and occurs primarily in Rp, 
i.e. Rp/Rz 8LS t —^ tc, consistent with recent exper- 
imental observations |17] . The same behavior occurs in 
lower- hybrid collapse in turbulent plasmas [T^j 130]. In 
the case of s-wave scattering, both Rp and Rz collapse 
at the same rate, with i?p (x (1 — t/tcY/'^, analogous to 
electrostatic Zakharov collapse [TH]. Furthermore, the 
fact that the TF parabolic scaling solutions give such a 
good description of the collapse indicates that, for the 
parameters considered, collapse is primarily a global ef- 
fect and proceeds through a quadrupole collective mode. 
From the GPE simulations a collapse time is defined in 
terms of a sudden onset of loss (red triangle) for K^> Q. 
Numerical simulations indicate that in the limit — > 
this coincides with the time at which p = Q (green circle) . 
Hence, our results indicate that GPE simulations with- 
out three-body losses can be used to infer an upper limit 
for the collapse time. For pancake shaped geometries 
(7^1) the TF and GPE solutions exhibit significant 
oscillations in Az before collapse, as seen in Fig. 2(c). 

From Eq. Q, using the GA, an upper bound for the 
collapse can be calculated (black squares in Fig. 2). For 
7^ — 0.1 the upper bound, is consistent with GPE re- 
sults. As we move away from the regime where the GA 
is valid. Figs. 2(b,c), the upper bound for the collapse 
time becomes inconsistent with the GPE results. 

Figure 3 presents the collapse times as kg is switched 
from fcso to ks, for various geometries, as evaluated via 
Eq. ^ (solid black), Eqs. ^\ (black circles), and GPE 
simulations with (red squares) and without (green -I-) 
three-body losses. For kg < k^ [left of the vertical dashed 
lines], the final state solutions are collapsed states, with 
the time it takes to collapse increasing, as kg approaches 
kg, from below. In all of the regimes presented we find 
that the TF approximation provides a good estimate for 
the upper bound of the collapse time. In contrast, Eq. 
([4|, is only consistent with the GPE simulations in the 
limit of very weak interactions. Note the appearance 
of a sudden step in the TF and GPE collapse times in 
Fig. 3(f). When the dipolar interactions dominate, the 
collapse is highly anisotropic and complete collapse oc- 
curs first in Rp. However, for large and attractive s-wave 
interactions, the collapse can become almost isotropic 
and for pancake-shaped systems complete collapse can 
occur first in i?^. This step represents the transition be- 
tween a complete collapse in Rz (left of the step) and Rp 
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FIG. 3: Collapse time according to Eq. Q (solid black curve), Eqs. (8j9 1 (black circles) and numerical integration of the GPE 
with (red squares) and without (green -I-, p=0) 3-body loss. For numerical solutions of the GPE with (without) loss we have 
used K3 = 2x 10"^* cmVs [E] {K3 = 0). (a-c) [(d-f)] fc^o = 60 [200] and kao = 30 [50], with (a) [(d)] 7^ = 0.1, (b) [(e)] 1.0 
and (c) [(f)] 10.0. The vertical dashed lines correspond to as derived from static solutions to the GPE. 



(right of the step). 

In summary we have presented two simplified mod- 
els of the collapse time of a DBEC, and compared them 
with exact numerical integration of dipolar GPE. When 
the atomic interactions are weak or attractive a GA for 
the initial DBEC can be assumed with an upper bound 
for the collapse time derived through a highly simpli- 
fied equation of motion for the radius [Eq. (4)]. In the 
opposing regime where interactions dominate, TF equa- 
tions of motion for the radii provide a convenient and 
approachable method to derive the time for global col- 
lapse. The validity of these two regimes is determined by 
the strength of the interactions and the aspect ratio of the 
trap. When dipolar interactions dominate, the collapse 
occurs in the plane perpendicular to the axis of polariza- 



tion. The excellent agreement, for the parameters con- 
sidered, between the TF dynamics and the full numerical 
simulations indicates that the collapse primarily occurs 
in a global manner and proceeds through a quadrupolar 
collective motion, consistent with the recent experimen- 
tal observations [111 HZ]. Finally, for oblate geometries, 
we observe two prominent deviations from this general 
behavior. Firstly, significant oscillations in can occur 
during collapse and secondly, for strong attractive s-wave 
interactions the collapse predominantly occurs in z rather 
than p. 
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